
Bienvenid@s a la segunda tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la segunda parte del curso, los cuales se enfocan principalmente en inferencia estadística, diseño de experimentos y test de hipótesis. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.
La tarea consta de una parte teórica que busca evaluar conceptos vistos en clases. Seguido por una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.
Slides de las clases:
Enlaces a videos de las clases:
Documentación:
En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.
Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 2 son las siguientes:
# Manipulación de estructuras
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
# Para realizar plots
library(ggplot2)
library(plotly)
##
## Adjuntando el paquete: 'plotly'
##
## The following object is masked from 'package:ggplot2':
##
## last_plot
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following object is masked from 'package:graphics':
##
## layout
# Manipulación de varios plots en una imagen.
library(gridExtra)
##
## Adjuntando el paquete: 'gridExtra'
##
## The following object is masked from 'package:dplyr':
##
## combine
Elimine eval=FALSE para ejecutar las celdas.
# Función para generar los estimadores
calcular_estimador <- function(n, sigma) {
# Simulación de las X_i ~ Bernoulli(0.5)
X <- rbinom(n, 1, 0.5)
# Ruido epsilon_sigma ~ N(0, sigma)
epsilon_sigma <- rnorm(1, mean = 0, sd = sigma)
# Estimador
p_gorro <- epsilon_sigma + (1/n) * sum(X)
return(p_gorro)
}
# sigma = 0
cte <- rep(0.5,1000)
x <- 1:1000
# Dataframe para almacenar los resultados
resultados_0 <- data.frame(n = numeric(), sigma = numeric(), p_gorro = numeric())
for (n in x) {
# Calcular el estimador usando la función
p_gorro <- calcular_estimador(n, 0)
# Almacenar los resultados
resultados_0 <- rbind(resultados_0, data.frame(n = n, sigma = 0, p_gorro = p_gorro))
}
# Graficar
ggplot(resultados_0, aes(x = n, y = p_gorro)) +
geom_line(color = "blue") +
geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
labs(title = "Estimador de p (sigma = 0)",
x = "n",
y = expression(hat(p)[n, 0])) +
theme_minimal()
# Sesgo para sigma = 0
sesgo_estimador_sigma_0 <- mean(resultados_0$p_gorro) - 0.5
print(sesgo_estimador_sigma_0)
## [1] -0.0006551944
# sigma = 1
cte <- rep(0.5,1000)
x <- 1:1000
# Dataframe para almacenar los resultados
resultados_1 <- data.frame(n = numeric(), sigma = numeric(), p_gorro = numeric())
for (n in x) {
# Calcular el estimador usando la función
p_gorro <- calcular_estimador(n, 1)
# Almacenar los resultados
resultados_1 <- rbind(resultados_1, data.frame(n = n, sigma = 1, p_gorro = p_gorro))
}
# Graficar
ggplot(resultados_1, aes(x = n, y = p_gorro)) +
geom_line(color = "green") +
geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
labs(title = "Estimador de p (sigma = 1)",
x = "n",
y = expression(hat(p)[n, 1])) +
theme_minimal()
# Sesgo para sigma = 1
sesgo_estimador_sigma_1 <- mean(resultados_1$p_gorro) - 0.5
print(sesgo_estimador_sigma_1)
## [1] 0.02413056
# sigma = 2
cte <- rep(0.5,1000)
x <- 1:1000
# Dataframe para almacenar los resultados
resultados_2 <- data.frame(n = numeric(), sigma = numeric(), p_gorro = numeric())
for (n in x) {
# Calcular el estimador usando la función
p_gorro <- calcular_estimador(n, 2)
# Almacenar los resultados
resultados_2 <- rbind(resultados_2, data.frame(n = n, sigma = 2, p_gorro = p_gorro))
}
# Graficar
ggplot(resultados_2, aes(x = n, y = p_gorro)) +
geom_line(color = "red") +
geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
labs(title = "Estimador de p (sigma = 2)",
x = "n",
y = expression(hat(p)[n, 1])) +
theme_minimal()
# Sesgo para sigma = 2
sesgo_estimador_sigma_2 <- mean(resultados_2$p_gorro) - 0.5
print(sesgo_estimador_sigma_2)
## [1] -0.006315325
# sigma = 4
cte <- rep(0.5,1000)
x <- 1:1000
# Dataframe para almacenar los resultados
resultados_4 <- data.frame(n = numeric(), sigma = numeric(), p_gorro = numeric())
for (n in x) {
# Calcular el estimador usando la función
p_gorro <- calcular_estimador(n, 4)
# Almacenar los resultados
resultados_4 <- rbind(resultados_4, data.frame(n = n, sigma = 4, p_gorro = p_gorro))
}
# Graficar
ggplot(resultados_4, aes(x = n, y = p_gorro)) +
geom_line(color = "purple") +
geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
labs(title = "Estimador de p (sigma = 4)",
x = "n",
y = expression(hat(p)[n, 1])) +
theme_minimal()
# Sesgo para sigma = 4
sesgo_estimador_sigma_4 <- mean(resultados_4$p_gorro) - 0.5
print(sesgo_estimador_sigma_4)
## [1] -0.1766147
Respuesta
En las preguntas 2, 3 y 4 deberá trabajar con el dataset
diabetes_prediction_dataset.
datos <- read.csv("diabetes_prediction_dataset.csv", header = TRUE)
head(datos)
Una forma de estimar la distribución de la media de una población es
a través de la realización de la sampling distribution de
una distribución cualquiera. El valor obtenido en cada sampling
distribution nos entrega un estimador de la media que posee una
determinada distribución de la población. Sabiendo esto, calcule la
sampling distribution de las columnas bmi y
blood_glucose_level, obteniendo el intervalo de confianza
de \(95\%\) para cada una de las medias
obtenidas desde la distribución. Para realizar este experimento
considere los siguientes puntos:
Hints:
bmi y blood_glucose_level a una misma
escala.sampling distribution podría serle
útil el comando sample().Respuesta:
# Definimos tamaños de muestreo
tamano_muestra = 100
n_muestras = 5000
# Inicializamos estructuras necesarias
list_mean_bmi = vector()
list_interval_bmi = vector()
list_typeCI_bmi = vector()
list_prop_bmi = vector()
sucess_bmi = 0
list_mean_bgl = vector()
list_interval_bgl = vector()
list_typeCI_bgl = vector()
list_prop_bgl = vector()
sucess_bgl = 0
#función para calcular desviación estándar
sd.p <- function(x){sd(x)*sqrt((length(x)-1)/length(x))}
# Obtenemos la media y desviación estándar de cada columna
mu_bmi <- mean(datos$bmi)
sd_bmi <- sd.p(datos$bmi)
mu_bmi
## [1] 27.32077
sd_bmi
## [1] 6.63675
mu_glucose <- mean(datos$blood_glucose_level)
sd_glucose <- sd.p(datos$blood_glucose_level)
mu_glucose
## [1] 138.0581
sd_glucose
## [1] 40.70793
alpha <- 0.05
# Sampling distribution, calculo del intervalo de confianza y proporción.
for (i in 1:n_muestras) {
### Para bmi ###
samples_bmi <- sample(datos$bmi, tamano_muestra, replace = TRUE)
sample_mu_bmi <- mean(samples_bmi)
sample_sd_bmi <- sd(samples_bmi)
# Calculo de intervalo de confianza
se_bmi <- sample_sd_bmi/sqrt(tamano_muestra)
error_bmi <- qnorm(1-alpha/2) * se_bmi
left_bmi <- sample_mu_bmi - error_bmi
right_bmi <- sample_mu_bmi + error_bmi
# Se guardan los resultados
list_mean_bmi[i] <- sample_mu_bmi
list_interval_bmi <- c(list_interval_bmi, left_bmi, right_bmi)
list_prop_bmi[i] <- sucess_bmi/i
# Se verifica si la media de la población se encuentra dentro del intérvalo de confianza
if (left_bmi <= mu_bmi && right_bmi >= mu_bmi) {
sucess_bmi <- sucess_bmi + 1
}
### Para blood_glucose_level ###
samples_glucose <- sample(datos$blood_glucose_level, tamano_muestra, replace = TRUE)
sample_mu_glucose <- mean(samples_glucose)
sample_sd_glucose <- sd(samples_glucose)
# Calculo de intervalo de confianza
se_glucose <- sample_sd_glucose / sqrt(tamano_muestra)
error_glucose <- qnorm(1 - alpha / 2) * se_glucose
left_glucose <- sample_mu_glucose - error_glucose
right_glucose <- sample_mu_glucose + error_glucose
# Se guardan los resultados
list_mean_bgl[i] <- sample_mu_glucose
list_interval_bgl <- c(list_interval_bgl, left_glucose, right_glucose)
list_prop_bgl[i] <- sucess_bgl / i
# Se verifica si la media de la población se encuentra dentro del intérvalo de confianza
if (left_glucose <= mu_glucose && right_glucose >= mu_glucose) {
sucess_bgl <- sucess_bgl + 1
}
}
# Generamos dataframes para plotear mas facilmente los datos.
# Convertimos los vectores de intervalos en matrices de dos columnas
list_interval_bmi_matrix <- matrix(list_interval_bmi, ncol = 2, byrow = TRUE)
list_interval_bgl_matrix <- matrix(list_interval_bgl, ncol = 2, byrow = TRUE)
# Generamos dataframes para graficar más fácilmente los datos
df_bmi <- data.frame(media = list_mean_bmi, ci_low = list_interval_bmi_matrix[, 1], ci_high = list_interval_bmi_matrix[, 2], prop_ci = list_prop_bmi)
df_bgl <- data.frame(media = list_mean_bgl, ci_low = list_interval_bgl_matrix[, 1], ci_high = list_interval_bgl_matrix[, 2], prop_ci = list_prop_bgl)
# Verificar si el intervalo de confianza contiene la media de la población
df_bmi$contiene_media <- (df_bmi$ci_low <= mu_bmi & df_bmi$ci_high >= mu_bmi)
df_bgl$contiene_media <- (df_bgl$ci_low <= mu_glucose & df_bgl$ci_high >= mu_glucose)
sucess_bmi
## [1] 4767
sucess_bgl
## [1] 4753
# Plot de Intervalos de confianza
# Graficar una muestra aleatoria de 100 medias e intervalos de confianza para BMI
set.seed(123) # Para reproducibilidad
sample_indices_bmi <- sample(1:n_muestras, 100) # Seleccionar 100 índices aleatorios
sample_bmi <- df_bmi[sample_indices_bmi, ]
# Crear gráfico para la muestra de BMI
fp_bmi <- ggplot(data = sample_bmi, aes(x = factor(1:100), y = media, ymin = ci_low, ymax = ci_high)) +
geom_pointrange() +
geom_hline(yintercept = mu_bmi, lty = 2) + # Línea horizontal para la media poblacional
coord_flip() + # Girar coordenadas
xlab("Muestra Aleatoria de Medias BMI") + ylab("Media (95% CI)") +
theme_bw() +
theme(axis.text.y = element_blank(), axis.ticks.y = element_blank()) # Opcional: limpiar etiquetas del eje y
print(fp_bmi)
# Plot de Intervalos de confianza
# Graficar una muestra aleatoria de 100 medias e intervalos de confianza para Glucose
set.seed(123) # Para reproducibilidad
sample_indices_bgl <- sample(1:n_muestras, 100) # Seleccionar 100 índices aleatorios
sample_bgl <- df_bgl[sample_indices_bgl, ]
# Crear gráfico para la muestra de Glucose
fp_bgl <- ggplot(data = sample_bgl, aes(x = factor(1:100), y = media, ymin = ci_low, ymax = ci_high)) +
geom_pointrange() +
geom_hline(yintercept = mu_glucose, lty = 2) + # Línea horizontal para la media poblacional
coord_flip() + # Girar coordenadas
xlab("Muestra Aleatoria de Medias Glucose") + ylab("Media (95% CI)") +
theme_bw() +
theme(axis.text.y = element_blank(), axis.ticks.y = element_blank()) # Opcional: limpiar etiquetas del eje y
print(fp_bgl)
# Para comparar estimaciones de bmi y blood_glucose_level
# Función que normaliza las variables a una escala entre 0 y 1
normalizar <- function(x) {
return((x - min(x)) / (max(x) - min(x)))
}
# Normalizamos las medias obtenidas de las muestras
list_mean_bmi_normalized <- normalizar(list_mean_bmi)
list_mean_bgl_normalized <- normalizar(list_mean_bgl)
# Creamos un data frame para facilitar el plotting
df_comparacion <- data.frame(
index = 1:n_muestras,
bmi = list_mean_bmi_normalized,
blood_glucose_level = list_mean_bgl_normalized
)
# Convertimos el data frame a formato largo sin reshape2
df_long <- data.frame(
index = rep(df_comparacion$index, 2),
Variable = rep(c("BMI", "Blood Glucose Level"), each = n_muestras),
Valor = c(df_comparacion$bmi, df_comparacion$blood_glucose_level)
)
# Graficamos usando ggplot2
library(ggplot2)
ggplot(df_long, aes(x = index, y = Valor, color = Variable)) +
geom_line() +
labs(title = "Comparación de Estimaciones Normalizadas de BMI y Glucosa en Sangre",
x = "Índice de Muestra",
y = "Estimación Normalizada") +
theme_minimal() +
theme(legend.title = element_blank())
# Plot de proporción de CI
p_prop_bmi <- ggplot(df_bmi, aes(x = 1:n_muestras, y = prop_ci)) +
geom_line(color = "blue") +
ggtitle("Proporción de Intervalos de Confianza que Contienen la Media - BMI") +
xlab("Iteraciones") + ylab("Proporción Acumulada") +
ylim(0, 1)
p_prop_bgl <- ggplot(df_bgl, aes(x = 1:n_muestras, y = prop_ci)) +
geom_line(color = "green") +
ggtitle("Proporción de Intervalos de Confianza que Contienen la Media - Glucose") +
xlab("Iteraciones") + ylab("Proporción Acumulada") +
ylim(0, 1)
# Mostrar las gráficas
print(p_prop_bmi)
print(p_prop_bgl)
Respuesta:
Es posible observar que la estimación para
blood_glucose_levelparece tener mejores resultados ya que a simple vista se observan más valores cercanos a la media real, en cambio parabmise observan valores mas dispersos y no tan cercanos a la media real por lo que se espera que este estimador tenga un mayor error. Sin embargo, lo observado puede deberse a otras razones como que los datos deblood_glucose_leveltienen una mayor dispersión de los datos y que por está razón se observan intérvalos de confianza más altos y por tanto pueda parecer que están más cercanos a la media real. Además es posible notar además al normalizar los estimadores deblood_glucose_levelybmique los estimadores de blood_glucose_level son mayores que parabmi, lo cual puede deberse a la naturaleza de los datos ya queblood_glucose_levelabarca datos en un mayor rango quebmi.
El objetivo de esta parte será obtener los parámetros de las
distribuciones de la media y la mediana de bmi. Para
responder la pregunta realice los siguientes puntos:
Cabe señalar que el método de máxima verosimilitud deberá ser programado por usted y no podrá utilizar librerías que entreguen el valor directo (por ejemplo la librería MASS).
Respuesta
# Obtenemos métricas de muestreo con repetición
set.seed(99)
means <- c()
medis <- c()
n_muestras <- 10000
tamano_muestra <- 100
for (i in 1:n_muestras) {
muestra <- sample(datos$bmi, tamano_muestra, replace = TRUE)
means <- c(means, mean(muestra))
medis <- c(medis, median(muestra))
}
# dataframe con las medias y medianas de las muestras
resultados <- data.frame(media = means, mediana = medis)
# Histograma media
ggplot(resultados, aes(x = media)) +
geom_histogram(binwidth = 0.1, fill = "blue", color = "white") +
ggtitle("Histograma de Medias Simuladas") +
theme_minimal()
# Histograma mediana
ggplot(resultados, aes(x = mediana)) +
geom_histogram(binwidth = 0.1, fill = "green", color = "white") +
ggtitle("Histograma de Medianas Simuladas") +
theme_minimal()
# Media
# función log likelihood
ll_plot <- function(mu, sigma) {
likelihood <- -sum(dnorm(means, mean = mu, sd = sigma, log = TRUE))
return(likelihood)
}
ll_plot <- Vectorize(ll_plot)
# Definir el espacio donde se va a evaluar ll_plot
mu_m <- seq(26, 29, length.out = 10) # Rango de valores para la media
sigma_m <- seq(0.01, 1, length.out = 10) # Rango de valores para la desviación estándar
likelihood_m <- outer(X = mu_m, Y = sigma_m, ll_plot)
filled.contour(x = mu_m, y = sigma_m, z = likelihood_m, xlab = expression(mu),
ylab = expression(sigma), color = topo.colors,
main = "Log Likelihood para la Media")
# Mediana
ll_plot_mediana <- function(mu, sigma) {
likelihood <- -sum(dnorm(medis, mean = mu, sd = sigma, log = TRUE))
return(likelihood)
}
ll_plot_mediana <- Vectorize(ll_plot_mediana)
mu_d <- seq(27, 28, length.out = 10)
sigma_d <- seq(0.01, 1, length.out = 10)
likelihood_d <- outer(X = mu_d, Y = sigma_d, ll_plot_mediana)
filled.contour(x = mu_d, y = sigma_d, z = likelihood_d, xlab = expression(mu),
ylab = expression(sigma), color = topo.colors,
main = "Log Likelihood para la Mediana")
likelihood_mean <- function(param) {
# Definimos los parametros de entrada de la funcion
mu <- param[1]
sigma <- param[2]
# Definimos la likelihood como la suma logarítmica de la función de densidad
log_likelihood <- -sum(dnorm(means, mean = mu, sd = sigma, log = TRUE))
return(log_likelihood)
}
# Agrrgue el rango donde operará nlminb
result_mean <- nlminb(objective = likelihood_mean,
start = c(mu = 27, sigma = 1), # Valores iniciales
lower = c(mu = 26, sigma = 0), # Límite inferior
upper = c(mu = 28, sigma = 2)) # Límite superior
# Mostrar los resultados
cat("Estimación por Máxima Verosimilitud para la media:\n")
## Estimación por Máxima Verosimilitud para la media:
cat("Mu óptimo:", result_mean$par[1], "\n")
## Mu óptimo: 27.32434
cat("Sigma óptimo:", result_mean$par[2], "\n")
## Sigma óptimo: 0.6633704
likelihood_med <- function(param) {
# Definimos los parametros de entrada de la funcion
mu <- param[1]
sigma <- param[2]
# Definimos la likelihood como la suma logarítmica de la función de densidad
log_likelihood <- -sum(dnorm(medis, mean = mu, sd = sigma, log = TRUE))
return(log_likelihood)
}
# Agrrgue el rango donde operará nlminb
result_med <- nlminb(objective = likelihood_med,
start = c(mu = 27, sigma = 1), # Valores iniciales
lower = c(mu = 26, sigma = 0), # Límite inferior
upper = c(mu = 28, sigma = 2)) # Límite superior
cat("Estimación por Máxima Verosimilitud para la mediana:\n")
## Estimación por Máxima Verosimilitud para la mediana:
cat("Mu óptimo:", result_med$par[1], "\n")
## Mu óptimo: 27.30507
cat("Sigma óptimo:", result_med$par[2], "\n")
## Sigma óptimo: 0.1025047
# Muestra de medias
mu_media <- 27.32434
sigma_media <- 0.663372
n_muestra <- 10000
muestras_media <- rnorm(n_muestra, mean = mu_media, sd = sigma_media)
df_muestras_media <- data.frame(muestras_media)
# Histogrmas
hist_media_simulada <- ggplot(df_muestras_media, aes(x = muestras_media)) +
geom_histogram(binwidth = 0.1, fill = "blue", color = "white") +
ggtitle("Histograma de Medias Simuladas") +
theme_minimal()
hist_media_simulada
# Muestra de medianas
mu_mediana <- 27.30507
sigma_mediana <- 0.1025046
n_muestra <- 10000
muestras_mediana <- rnorm(n_muestra, mean = mu_mediana, sd = sigma_mediana)
df_muestras_mediana <- data.frame(muestras_mediana)
# Histogrmas
hist_mediana_simulada <- ggplot(df_muestras_mediana, aes(x = muestras_mediana)) +
geom_histogram(binwidth = 0.1, fill = "green", color = "white") +
ggtitle("Histograma de Medianas Simuladas") +
theme_minimal()
hist_mediana_simulada
Al analizar las distribuciones de medias y medianas obtenidas por máxima verosimilitud, se observa lo siguiente:
Medias: La estimación de máxima verosimilitud captura adecuadamente la variabilidad observada en los datos empíricos. Esto se refleja al comparar los histograma de las medias empíricas y el de las medias simuladas, lo que indica que la distribución normal estimada es un buen ajuste para los datos.
Medianas: En el caso de las medianas, no se observa el mismo comportamiento. El histograma de las medianas empíricas muestra una fuerte concentración en un solo valor, lo que indica una baja variabilidad en los datos originales. Esto difiere con las medianas simuladas, que presentan mayor variabilidad, propia de una distribución más normal.
La diferencia en la variabilidad entre las medianas empíricas y las simuladas sugiere que la distribución de los datos empíricos tiene características que no son bien capturadas por la estimación normal, posiblemente debido a una concentración excesiva de valores en un punto específico.
En esta pregunta se error e inetrvalo de confianza para la varianza
de la columna HbA1c_level.
Las actividades por realizar son las siguientes:
Respuesta:
# Bootstrap
B <- 5000
largo <- length(datos$HbA1c_level)
output <- c()
for (it in 1:B) {
sample_data <- sample(datos$HbA1c_level, largo, replace = TRUE)
output[it] <- var(sample_data)
}
# Valor real de la varianza
real_variance <- var(datos$HbA1c_level)
real_variance
## [1] 1.146339
# Gráfico de puntos
plot(output, main = "Varianzas obtenidas con Bootstrap",
ylab = "Varianza", xlab = "Muestra de Bootstrap")
abline(h = real_variance, col = "red", lty = 2) # Línea roja para la varianza real
legend("topright", legend = "Varianza real", col = "red", lty = 2)
# Histogrma
hist(output, main = "Histograma de varianzas de Bootstrap",
xlab = "Varianza", breaks = 30, col = "lightblue")
abline(v = real_variance, col = "red", lty = 2) # Línea roja para la varianza real
legend("topright", legend = "Varianza real", col = "red", lty = 2)
# Obtenemos error e intervalos de confianze
se_vars <- sd(output)
z_a2 <- qnorm(1 - alpha/2)
alpha <- 0.05
# límite inferios
l.CI <- mean(output) - z_a2 * se_vars
# límite superior
u.CI <- mean(output) + z_a2 * se_vars
sprintf('El intervalo de confidencia de 95%% de las varianzas es: (%.3f,%.3f)', l.CI , u.CI)
## [1] "El intervalo de confidencia de 95% de las varianzas es: (1.136,1.157)"
sprintf('El SE de la varianzas es: (%.3f)', se_vars)
## [1] "El SE de la varianzas es: (0.005)"
Respuesta:
Para el gráfico de puntos para la varianza obtenida con boostrap se puede ver que los puntos están mayormente concentrados en la línea punteada roja, la cual corresponde a la varianza real. Además se nota que hay puntos outliers presentes pero en menor cantidad.
Para el histrograma de la varianza se cómo claramante las varianzas obtenidas con Bootstrap tienen una distribución cercana a la normal centrada en el valor de la varianza real.
Se puede inferir que al aplicar el método bootstrap para error estándar de la varianza, estos resultados tendrán una distribución cercana a la normal centrada en el valor de la varianza real, más aún en el histograma ya se puede apreciar cuál será un valor aproximado de el 95% de intervalo de confianza de la estimación.
El siguiente código genera una regresión lineal de bmi
con respecto a age, es decir, una función lineal de
age, \(y(age) = b + m\cdot
age\), que pretende determinar el valor de bmi.
linearMod <- lm(bmi ~ age, data=datos)
print(linearMod)
##
## Call:
## lm(formula = bmi ~ age, data = datos)
##
## Coefficients:
## (Intercept) age
## 23.15536 0.09945
m <- linearMod$coefficients["age"]
b <- linearMod$coefficients["(Intercept)"]
ggplot() +
geom_point(aes(x=datos$age, y=datos$bmi)) +
geom_line(aes(x=datos$age, y=datos$age*m+b), color="red") +
ggtitle("Regresión lineal") +
ylab("bmi") +
xlab("Age")
Repita el proceso de la pregunta 4 para estimar el error e intervalos de confianza de los parámetros de la regresión (\(m\) y \(b\)). ¿Qué indican los resultados de Bootstrap?¿Un valor bajo en el error significa que la regresión es buena?
# Bootstrap
B <- 500
largo <- 1000
output_m <- numeric(B)
output_b <- numeric(B)
for (it in 1:B) {
# Tomar una muestra con reemplazo de los datos
sample_indices <- sample(1:largo, largo, replace = TRUE)
sample_data <- datos[sample_indices, ]
# Ajustar la regresión lineal con la muestra
bootstrap_model <- lm(bmi ~ age, data = sample_data)
# Guardar los coeficientes en las salidas correspondientes
output_m[it] <- bootstrap_model$coefficients["age"]
output_b[it] <- bootstrap_model$coefficients["(Intercept)"]
}
# Mostrar algunos resultados del bootstrap
head(output_m)
## [1] 0.10276581 0.09021273 0.09894292 0.09804121 0.10033194 0.09843114
head(output_b)
## [1] 22.70921 22.90499 22.90196 23.14547 22.67139 22.92839
# Obtenemos error e intervalos de confianze
se_m <- sd(output_m)
z_a2 <- qnorm(1 - alpha / 2)
alpha <- 0.05
# límite inferios
l.CI <- mean(output_m) - z_a2 * se_m
# límite superior
u.CI <- mean(output_m) + z_a2 * se_m
sprintf('El intervalo de confidencia de 95%% de las varianzas es: (%.3f,%.3f)', l.CI , u.CI)
## [1] "El intervalo de confidencia de 95% de las varianzas es: (0.080,0.113)"
sprintf('El SE de la varianzas es: (%.3f)', se_m)
## [1] "El SE de la varianzas es: (0.008)"
# Obtenemos error e intervalos de confianze
se_b <- sd(output_b)
z_a2 <- qnorm(1 - alpha / 2)
alpha <- 0.05
# límite inferios
l.CI <- mean(output_b) - z_a2 * se_b
# límite superior
u.CI <- mean(output_b) + z_a2 * se_b
sprintf('El intervalo de confidencia de 95%% de los sesgos es: (%.3f,%.3f)', l.CI , u.CI)
## [1] "El intervalo de confidencia de 95% de los sesgos es: (22.208,23.759)"
sprintf('El SE del sesgo es: (%.3f)', se_b)
## [1] "El SE del sesgo es: (0.396)"
Respuesta
Error estándar bajo: Indica que los coeficientes son consistentes a través de las muestras, lo que sugiere un modelo de regresión estable.
Un error bajo no garantiza que el modelo sea bueno en términos predictivos, pero lo hace fiable de que observa relaciones en los datos actuales.
A work by CC6104